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ABSTRACT 

Context. We interpret solar flares as events originating from active regions that have reached the Self Organized Critical state, by 
using a refined Cellular Automaton model with initial conditions derived from observations. 

Aims. We investigate whether the system, with its imposed physical elements, reaches a Self Organized Critical state and whether 
well-known statistical properties of flares, such as scaling laws observed in the distribution functions of characteristic parameters, are 
reproduced after this state has been reached. 

Methods. To investigate whether the distribution functions of total energy, peak energy and event duration follow the expected scal- 
ing laws, we first apply a nonlinear force-free extrapolation which reconstructs the three-dimensional magnetic fields from two- 
dimensional vector magnetograms. We then locate magnetic discontinuities exceeding a threshold in the Laplacian of the magnetic 
field. These discontinuities are relaxed in local diffusion events, implemented in the form of Cellular Automaton evolution rules. 
Subsequent loading and relaxation steps lead the system to Self Organized Criticality, after which the statistical properties of the 
simulated events are examined. Physical requirements, such as the divergence-free condition for the magnetic field vector, are approx- 
imately imposed on all elements of the model. 

Results. Our results show that Self Organized Criticality is indeed reached when applying specific loading and relaxation rules. Power 
law indices obtained from the distribution functions of the modeled flaring events are in good agreement with observations. Single 
power laws (peak and total flare energy) as well as power laws with exponential cutoff and double power laws (flare duration) are 
obtained. The results are also compared with observational X-ray data from GOES satellite for our active-region sample. 
Conclusions. We conclude that well-known statistical properties of flares are reproduced after the system has reached Self Organized 
Criticality. A significant enhancement of our refined Cellular Automaton model is that it commences the simulation from observed 
vector magnetograms, thus facilitating energy calculation in physical units. The model described in this study remains consistent with 
fundamental physical requirements, and imposes physically meaningful driving and redistribution rules. 

Key words. Sun: corona - Sun: flares - Sun: photosphere 
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C> 1 . Introduction rence, based on the assumption that the solar corona is in a sta- 

y ~ ' ' tistically stable Self Organized Critical (SOC) state. In this con- 

t-H ; Solar flares are transient energy release events above so- text> ARs are per ceived as nonlinear dissipative dynamical sys- 



X 



lar Active Regions (ARs). Populations of flares are known tems? externally driven by the photospheric velocity field. Due 
to exhibit robust statistical properties, which have been re- t0 random shuffling of the coronal loops' footpoints in the pho- 



peatedly identified in numerous observations. In partial- tosphere, localized instabilities are generated, which are respon 

U • lar, specific flare parameters have been consistently found sible for the fragme nted energy release in the solar corona. The 

W.-to follow robust power laws with indices lying in well- magnetic energy release simulated via Cellular Automaton (CA) 

define d ranges. More specifically, a series of flare observa- mod eling led to avalanche-like events. This model allows insta- 

tions (IDatlowe et al. 19741 I Di L et_aLJ.9 84 l Sturrock et al. 19841 bilitieSj s i mu i ating curren t sheet disruption and Ohmic dissipa- 

IDennis 19851 IVilmer 19871 |Crosby et al. 1993] | Biesecker 1994] tion> when a certain current density thres hold is exceeded. An 



Bromundetal.T29Sl |Polygiannakis et al. 2002| l report that the enhancement of the original SOC concept with respect to the 



distribution functions of peak flux, total energy and event du- instab ii ity criteria and the co rres ponding relaxation was intro- 

ration exhibit well-formed scaling laws with exponents in the duced by [Viihos et al. (1995)| and |Ueorgoulis et al. (1995)[ Both 

ranges of (-1.59,-1.80), (-1.39, -1.50) and (-2.25, -2.80) re- SU gg es ted that the initial instability may trigger secondary ones, 

spectively. mus a ff e cting sites beyond the closest vicinity of the original 

This consistency of the power law indices identified in nu- event. Non-lo cal treatment between flaring eleme nts was also 

merous independent studies stimulated a new phenomenological attempted by |MacKinnon & Macpherson (1997)] In addition, 



approach in reprod ucing and modeling the stati stical behaviour |Georgoulis & Vlahos (1996)| constructed a refined Statistical 



of flaring activity. |Lu & Hamilton (1991 )| and |Lu et al. (1993)1 Flare Model, including both isotropic and anisotropic relaxation 
were the first to construct a simple model of solar flare occur- 
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mechanisms as well as extended instability criteria. This kind 
of modeling produced a double power-law scaling behaviour: 
the flatter power law resembled intermediate and large flares, 
whereas the steeper one described low-energy events. An addi- 
tional enhancement of this model lay in the external driver simu- 
lation. In the mentioned study, the driver of the system followed 
a power law itself, thus mimicking the instabilities triggered by 
the emerging magnetic flux from the convection zone in addition 
to photospheric shuffling. Further extensions were introduced by 
Georgoulis & Vlahos (1998), who presented a systematic study 
of the power law indices' variability as a function of the driver's 
properties. In this refined Statistic Flare Model, Georgoulis & 
Vlahos attempted to model the stresses which are built-up ran- 
domly within ARs through a highly variable, inhomogenous ex- 
ternal driver. Although clearly deviating from the initial SOC 
models, the robust scaling laws in the flares' distribution func- 
tions survived. Isliker et al. (2000. 2001 ) tried for the first time 
to associate the classical CA models' components with physical 
variables. Magnetohydrodynamic (MHD) and CA approaches 
were connected through the physical interpretation of numerous 
CA elements, like the grid-variable, the time step, the spatial dis- 
creteness, the energy release process and the role of diffusivity. 
This study revealed several inconsistencies of the CA modeling, 
such as the uncontrolled value of the magnetic field divergence 
(V • B) and the non-availability of secondary variables, such as 
the current density and the electric field. Such weaknesses were 
treated by the Extended CA model (X-CA) introduced by Isliker 
et al. 1(2002)| 

In this study we present a model which adopts the Lu & 
Hamilton (1991) approach as starting point and to which sev- 
eral enhancements are made towards a more physical CA model 
that integrates various aspects of observed ARs and flares: first 
and foremost, the initial boundary and initial conditions stem 
from observed vector magnetograms. This allows us to per- 
form calculations in physical units, in direct comparison with 
observations (see for example the respective restrictions pre- 
sented in Georgouli s et al. (2001)| >. Time remains the only quan- 
tity expressed in arbitrary model units, as the photospheric vec- 
tor magnetogram does not change during the simulation. An ad- 
ditional feature is that during the whole process (initial load- 
ing, relaxation of magnetic discontinuities and further driving) 
the requirement V • B ^ is explicitly imposed. For this pur- 
pose we have used a nonlinear force-free extrapolation method 
to generate the initial conditions from observed magnetograms 
and impose instability criteria related to actual physical pro- 
cesses. The magnetic field relaxation in the CA model follows 
the |Lu & Hamilton (1991)| principles. The driving process is also 
designed to obey specific rules which do not violate known phys- 
ical processes in the corona. The structure of this work is as fol- 
lows. Section 2 describes the data used in this study along with 
the necessary corrections imposed on them. Section 3 explains in 
detail all the modules comprising our model: first the extrapola- 
tion technique (EXTRA) along with the discontinuities' identifi- 
cation (DISCO) modules. Furthermore, the magnetic field relax- 
ation module (RELAX) and finally the driving module (LOAD) 
is presented, which may trigger further instabilities in the sim- 
ulated AR, following rules which mimic specific physical pro- 
cesses. Section 4 presents our results and discusses our findings. 
Finally, Sect. 5 summarizes our conclusions. 

2. Data-set 

Nonlinear force-free extrapolation techniques require vector 
magnetograms that are not as widely available as conventional 



line-of-sight magnetograms. Here we have created a database 
of 1 1 different AR vector magnetograms from the University of 
Hawaii Imaging Vector Magnetograph (IVM). 

IVM obtains Stokes images in photospheric lines with 7pm 
spectral resolution, l.larcsec spatial resolution (~ 0.55 arcsec 
per pixel) over a field of A.larcmin 2 and polarimetric preci- 



sion of 0.1% ( Mickey et al. 1996). We use both fully-inverted 
and quick-look IVM data. Quick-look data have been ob- 
tained from the IVM Survey Data archive (available online at 
http : //www . cora . nwra . com/ivm/IVM- SurveyData). The 
quick-look data reduction differs from the complete inversion in 
that it uses a simplified flat-fielding approach, takes no account 
of scattered or parasitic light, and no correction is attempted for 
seeing variations that occur during the data acquisition. 

In this study we use 1 fully inverted and 10 quick-look 
IVM vector magnetograms. To remove the intrinsic azimuthal 
ambiguity of 180° we use the Non-Potential magnetic Field 
Calculation (NPFC) method of Georgoulis (2005). For compu- 
tational convenience we further rebin the disambiguated magne- 
tograms into a 32 x 32 regular grid. 

3. The model 

Our model consists of 4 separate modules. First we apply the 
Wiegelmann (2008) optimization algorithm to our vector mag- 
netograms in order to nonlinearly extrapolate the magnetic field 
from the photospheric boundary (module "EXTRA"). We thus 
construct a three-dimensional (3d) 32 x 32 x 32 cube, within 
which the magnetic field is determined. Second, we identify the 
sites within our cubic grid that exceed a threshold in the mag- 
netic field Laplacian (module "DISCO"). If unstable sites are 
found, we force the vicinity of the unstable location to undergo a 
magnetic-field restructuring. This redistribution is governed by 
specific rules, which do not violate basic physical laws. Under 
suitable conditions, the onset and relaxation of an initial instabil- 
ity may trigger a cascade of similar events in an avalanche-type 
manner. It is clear, therefore, that the wider vicinity, up to the en- 
tire system, may participate in this process. Module "RELAX" 
handles the field redistribution triggered by both the primary and 
subsequently triggered instabilities. The whole avalanche, com- 
prised of a seed, and subsequently triggered instabilities, is con- 
sidered as one single flare. After complete relaxation, we further 
drive the system via the "LOAD" module. There, a randomly- 
selected grid site receives a random magnetic field increment. 

3.1. "EXTRA": a nonlinear force-free extrapolation module 
The first step is the extrapolation of the photospheric mag- 



netic fields. As explained in |Dimitropoulou et al. (2009) 



physically meaningful treatment is the nonlinear force-free 
(NLFF) field extrapolation. Our method of choice is based 
on the optimization technique introduced by Wheatland 
et al. (2000) and further developed by T. Wiegelmann 



and collaborators (Wiegelmann 2004 Wiegelmann et al. 2006 



Wiegelmann 2008). This technique reconstructs force-free mag- 
netic fields from their boundary values by minimizing the 
Lorentz force and the divergence of the magnetic field vector 
in the extrapolation volume: 



= f w(x,y,z)[\B\- 2 \(VxB)xB\ 2 + \V -B\ 2 ]d 3 x 
Jv 



(1) 



In the above functional, w(x, y, z) is a weighting function and 
V denotes the extrapolation volume. A force-free state is reached 
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when L — > for w > 0. For w(x, y,z) = 1 the magnetic field 
must be available on all 6 boundaries of our cubic box for the 
optimization algorithm to work. However, photospheric vector 
magnetograms pertain only to the bottom boundary, whereas the 
magnetic field vector on the top and lateral boundaries is un- 
known. The weighting function is thus used in order to mini- 
mize the dependence of the interior solution from the unknown 
boundaries. In this study we introduce a buffer zone of 10 grid 
points expanding to the lateral and top boundaries of the compu- 
tational box. We then choose w(x, y, z) = 1 in the inner domain 
and let w drop to with a cosine-profile in the buffer zone to- 
wards the lateral and top boundaries of the computational box. 
This technique was first described by Wiegelmann (2004). 

An additional useful attribute of Wiegelmann's NLFF field 
extrapolation code is the preprocessing option it offers. As the 
photospheric magnetic field is in principle inconsistent with the 
force-free approximation, a preprocessing procedure was devel- 
oped by Wiegelmann et al. (2006) in order to drive photospheric 
fields closer to a NLFF field equilibrium. Preprocessing mini- 
mizes the forces and torques in the system, thus satisfying the 
force-free requirements more closely. 

Although NLFF extrapolation methods have been greatly 
improved over recent years, such models still include numerous 
uncertainties (DeR osa et al. 2009V Additional constraints stem 
from the measurements (signal-to-noise ratio, inadequate reso- 
lution of the 180° ambiguity) or from physical origins (variation 
in the line formation height, the non-force-free nature of the pho- 
tospheric vector magnetograms), which are not adequately han- 
dled in the course of the extrapolation. Such uncertainties are 
unavoidably conveyed to our simulations. 

3.2. "DISCO": a module to identify magnetic-field instabilities 

We assume that instabilities occur if the magnetic field stress 
exceeds a critical threshold. For every site r within our grid we 
calculate the magnetic field stress G m ,(r) as 

Gavto = \G ar (r)\ 

where 

G m {r) = B(f)-±- n Y Jnn B nn {r) 

In the above definitions nn is the number of nearest neigh- 
bors for each site r and B n „(r) is the magnetic field vector of 
these neighbors. Depending on the location of each site within 
the volume, the number of nearest neighbors nn can be nn = 
3,4,5,6. The physical reason for selecting this criterion lies in 
the fact that large magnetic stresses favor magnetic reconnec- 
tion in three dimensions, even in the absence of null points 
dPriest et al. 2003b . 

Mathematically, it can be shown that the selection of G av as 
the critical quantity in our model relates to the diffusive term 
of the induction equation (see |Isliker et al. (1998)| for a detailed 
discussion). Let us write the induction equation in the form: 

^ = V x (V x B) + riV 2 B (2) 
of 

where V is the plasma velocity and rj is the resistivity. The 
Laplacian of the magnetic field V 2 B(r) can be written as follows: 

V 2 B(r) = (V 2 B A )i + (V 2 B v )j + (V 2 fi,)k, 



where r — (z, j, k). Letting m = x,y,z we obtain 

gf " ~ ~&y2(Bm iJ+] j. + B mi Ht — 2B m . jt ) 

— ~^{B mijm + B mijkl — 2B,„ ijk ) 

adopting a central finite-difference scheme and using the 
general case of a grid point having 6 nearest neighbors (nn = 6). 
Further assuming Ax = Ay — Az = 1 (the grid-size) we have: 

V 2 B m (r) = *fjfa + *jf> + *fjfa , Zm Bmm - nnB m . 

which yields V 2 B(r) as follows: 

V 2 B(r)^ZnnBnn(r)-nnB(r) 

From the definition of the critical quantity G av {r) it follows 
that: 

V 2 B(r) - -nnG av (r) (3) 

Therefore the critical quantity G m ,(r) relates directly to the 
Laplacian V 2 B. The resistivity in the solar corona is almost zero 
everywhere except in regions where the discontinuities (and the 
local currents) reach a critical value. In these regions current- 
driven instabilities will enhance the resistivity by many orders 
of magnitude and the second term in equation (2) will become 
dominant. The convective term V x (V x B) of equation (2) will 
be further discussed in section 3.4, where the driving module 
"LOAD" is described. 

There are several ways to determine the threshold value for 
the critical quantity, above which a site is considered unstable: 

1 . We apply a histogram method, by constructing the histogram 
of the G av values in our grid. We then fit a Gaussian to 
this histogram and define the threshold G cr as the field 
stress value, above which the histogram deviates from the 
Gaussian. 

2. We define the threshold value (G cr ) for the whole grid, as 
the maximum G av , m „ value throughout our volume, slightly 
decreased: 

G cr - G aVmax {\ - s) 
where s « 1 

3. We define the threshold value (G cr (z)) per height z , as the 
maximum G aVmm value for each specific height, slightly de- 
creased: 

G cr (z) = G aVm Jz)(l - S) 

where s « 1 

4. We define the threshold value as a function of height z, e.g.: 

G cr {z) = G fllWr (l - s)exp(-z) 
where s « 1 

Here we present the results produced by the first (histogram) 
method, which yielded G cr = 1QG for our sample, and shortly 
refer to the other threshold alternatives in Sect.4. Every site 
r = (i,j,k) for which the inequality G av . > G cr is satisfied, is 
considered unstable and undergoes magnetic field restructuring 
under the rules implemented in the "RELAX" module. Notice 
that, given the definition of the critical threshold, instabilities 
sometimes occur even from the first iteration, after constructing 
the NLFF fields. This, however, does not incur any qualitative 
impact on the evolution of the system toward SOC, or the statis- 
tical results of the simulation. 
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3.3. "RELAX": a redistribution module for magnetic energy 



In case the instability criterion G aVijk > G cr is met for a 
specific site i,j, k, then the vicinity of the unstable location 
undergoes a field restructuring, which follows the rules of 
|Lu& Hamilton (1991)| 



B + {r) -» B(r) - -G av {r) 



B + nn (r) -» B m (r) + -G av {r) 



(4) 



(5) 



where the superscript + denotes the field components after 
the redistribution. 

At this point it is important to investigate whether the 
redistribution rules as defined here violate basic physical laws, 
such as the zero-divergence requirement for the magnetic field. 
The initial magnetic configuration satisfies approximately the 
condition V ■ B — 0, as the field has been reconstructed using an 
NLFF field extrapolation. The question is whether the magnetic 
field after the redistribution imposed by rules (4)-(5) still 
satisfies the same demand (V • B + = 0). Taking the divergence 
of B + (r) and its neighbors B* (r) we find from relations (4) and 
(5) respectively: 

V ■ B + (r) - V • B(r) - f V ■ G av (r) 

V ■ B + nn (r) - V ■ B m (r) + ±V ■ G av (r) 
From the definition of G av {r) we now have: 

V • G m (r) = V • B(r) - • B m (r) 
Substituting this to the above we find: 



V • B + (r) - iv ■ B(r) - J-V ■ B m (r) 

7 Inn 

V ■ B+ m (r) - iv • B(r) + -Lv ■ B m {r) 

I Inn 



(6) 



(7) 



As V ■ B(r) ^ V ■ B„„(r) ^ from our first iteration (ex- 
trapolated fields), we find V • B + (r) * V • B + m (r) 0. Thus, the 
redistribution of the magnetic field maintains the divergence-free 
condition. 

Isliker et al. (1998) showed that the redistribution rules 
(4) and (5) implement local diffusion and after redistribution, 
Go,,(r) st 0, so the instability at location r has been relaxed. 

3.4. "LOAD": the driver 

After the system is completely relaxed, we introduce a driving 
mechanism that adds a magnetic field increment 6B(r) at one 
randomly selected site r within our grid. The driving process 
complies with the following conditions: 

1. B{r)-SB{r) = (8) 
This condition implies that the magnetic field increment is 
always perpendicular to the existing magnetic field B(r) at 
the randomly selected site r. Figure 1 provides a sketch of the 
suggested situation, depicting the directions of the plasma 
velocity V, the magnetic field B and the perpendicular mag- 
netic field increment 6B. We note that the condition de- 
scribed by Eq. (8) is compatible with two physical scenar- 
ios: (a) that Alfven waves may have been excited locally, or 



magnetic 
field-line 




photosphere 



Fig. 1. Typical configuration of a magnetic loop anchored in the 
photosphere. The magnetic field vector B is perpendicular to an 
assumed plasma outflow velocity V. The model driver requires 
that the magnetic field increments 6B are always perpendicular 
to the existing magnetic field B. 



(b) that, according to the convective term V x (V x B) of the 
induction Eq. (2), a magnetized plasma upflow occurs in the 
AR, out from the photosphere. 

2. ^= e ,,<l ' (9) 

\B(r)\ 

This is a typical condition known to allow the system 
reach the SOC state, without the latter being influenced 
by the loading process (IBak et al. 19871 1. As also shown by 
|Lu & Hamilton (1991)| decreasing the driving rate by mak- 
ing the magnetic field increments even smaller, increases the 
average time between subsequent events. For the results pre- 
sented here we have used a fixed e = 0.3. 

3. V • (B(r) + 6B(r)) = (10) 
This condition should guarantee that the divergence of the 
magnetic field is approximately kept to zero during the 
loading process, as it was done during the redistribution of 
the magnetic field (RELAX module). In order to implement 
the condition, a first-order, left finite-difference scheme is 
used. In this way , however, condition (10) does not provide 
adequate guarantee for a divergence-free magnetic field in 
the selected site's vicinity. This is a known problem, which 
can be tackled by working with the vector potential A, with 
V x A — B, instead of the magnetic field B directly (see 
e.g. |Lu et al. (1993)| [Galsgaard (1996)| Isliker et al. (2000), 
Isliker et al. (2001)). Because our study uses observed vector 
magnetograms as initial conditions, we naturally work with 
the known magnetic fields, rather than the unknown vector 
potential. Thus, equation (10) only provides a low-order 
approximation towards a divergence-free magnetic field. 
To monitor how effective condition (10) is, we introduce a 
"Weighted Nabla Dot B" (WNDB) monitoring parameter, 
as follows: 



WNDB 



g-gl 



By definition, WNDB is a dimensionless quantity, lying in 
the range < WNDB < 1 . Monitoring WNDB during our 
simulation will provide evidence on whether condition (10) 
can be considered adequate for keeping the magnetic field 
within our volume approximately divergence-free. In the fol- 
lowing, we will tolerate a departure from zero of up to 20% 
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(WNDB < 0.2) for a still a roughly divergence-free magnetic 
field. 

3.5. Model parameters 

Linking the above-mentioned modules in one consistent simula- 
tion, we construct a relatively simple algorithm and we monitor 
the flare duration, the peak energy and the total energy, the dis- 
tribution functions of which we intend to compare against those 
of observational data. If an instability is identified (DISCO) - 
either in the original magnetic configuration generated by the 
initial extrapolated magnetogram (EXTRA) or due to an incre- 
ment SB randomly added at a grid point (LOAD) - the possi- 
ble chain of instabilities that follows is left to completely relax 
(RELAX) before an additional magnetic field increment is ran- 
domly placed (LOAD), possibly causing a new instability. This 
rule takes into account the observational fact that the lifetime of 
a flare is much smaller than the evolution timescale of an AR. 
Successive grid scans may be required for an instability to be 
completely relaxed. Each scanning corresponds to one timestep. 
Therefore, the relaxation of an event may be accomplished in 
more than one timesteps. Each loading event according to the 
equation set (8)-(10) commences a new iteration. The duration 
of an event is defined as the total number of timesteps the event 
lasted, from its onset until its complete relaxation. The accumu- 
lated released energy during the event provides the total energy 
of an event, whereas the peak during an event yields the peak 
energy/luminosity of an event. 

The simulation results presented in the next section have 
been performed using a 32x32x32 cubic grid with "open" 
boundaries in the relaxation events (see Isliker et al. (2001) for 
a detailed discussion on open boundary conditions). Each sim- 
ulation is driven for 3xl0 5 iterations, which equals the times 
that LOAD module is being called during the simulation. This 
mechanism allows the production of multiple subsequent flares 
in each AR. In all cases the critical threshold G C r 

was kept fixed 

and equal to 10G. 

4. Results 

Applying our flare simulation model to our 1 1 -event-database, 
we find that in all cases the simulated ARs reached the SOC 
state. An indication on whether and when the SOC state has 
been reached is obtained by monitoring the quantity G av , namely, 
the volume average of the critical quantity G av . During the con- 
tinuous driving of the system and the subsequently generated 
avalanches, G av increases gradually. When the system reaches 
the SOC state, G av stabilizes around a value which depends on 
the system's characteristics. For the loading method used in our 
model (new magnetic field increments are only added when a 
previously triggered avalanche has decayed), the value around 
which G av stabilizes is slightly lower than the threshold value 
G cr . A second indication that the system has reached the SOC 
state is that the total energy of the system tends to an asymp- 
totic value. This is because SOC is a statistically stationary state. 
Figure 2 shows the G av value over 3xl0 5 timesteps for AR 10570. 
G av is constantly increasing up to timestep 1.4xl0 5 , thereafter 
stabilizing at ~ 9.80G < G cr = IQG. Similarly, Fig. 3 shows the 
logarithm of the total magnetic energy throughout the volume 
Etotaft after each scan of the grid for possible redistributions. 
Following G av , E tota f, increases until a stable state is reached. 
The stabilization of both G av and E tota f t is a solid indication that 
SOC has been reached for AR 10570. The same behavior is seen 
for all ARs included in our sample. 
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Pi I 

1x10 s 2xl0 3 3x10 s 4x10 s 

Timesteps 

Fig. 2. Average Laplacian G~ av over the grid for 3x10 s timesteps 
for A/?10570. G av increases gradually until timestep 1.4xl0 5 , af- 
ter which the SOC state is reached. 

33.30C — ' — ' — ' — ' — i — ' — ' — ' — ' — i — ' — ' — ' — ■— ] 




o 

33.10 : 



33.00 h ■ 

5.0xl0 5 l.OxlO 6 1.5xl0 6 

# Redistributions 

Fig. 3. Diagram of log 10 (2? to<B / ( ) after each redistribution for 
AR 10570. Like G av , E tota f, increases gradually until a stable state 
is reached. 

SOC is generally characterized by intermittent transport 
events (avalanches), whose sizes range from very small (a sin- 
gle neighborhood) up to comparable to the system size. Power- 
law frequency distributions describe the parameters of these 
avalanches. It is thus reasonable to expect that since all 1 1 ARs 
in our sample have reached the SOC state under the imposed 
driving rules, they should all produce distribution functions for 
the flare duration, peak energy and total energy, which either 
follow pure power laws or functions including a power law part 
(e.g. power laws with exponential rollover). The functions tested 
against the model results for all flare parameters (flare dura- 
tion, peak energy, total energy) were single power laws, double 
power laws, power laws with exponential rollover, and exponen- 
tial functions. In order to define the best-fitting function per case, 
we made least square fits and performed chi-square goodness-of- 
fit tests. 

Figures 4 and 5 are typical examples of our general results. 
Figure 4 depicts the distribution functions of duration (Fig. 4a), 
peak energy (Fig. 4b) and total energy (Fig. 4c)) for A/?10050. 
The duration distribution follows a double power law with index 
-1.80 + 0.18 for the flatter part and -4.03 ± 0.29 for the steeper 
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Table 1. Power-law indices and respective chi-square probabilities derived by the best-fit functions for the flare duration for the 1 1 
ARs comprising our sample. 



FLARE DURATION (model) 

Double Power Law fit Power Law with 





Flat PL 


Steep PL 


Exponential Rollover 


AR 


PL Index 


Probability 


PL Index Probability 


PL Index 


Probability 


9415 








-1.42 ± 0.18 


0.95 


9635 


-2.29 ±0.19 


0.96 


-5.28 ± 0.42 0.95 






9661 








-0.26 ± 0.05 


0.94 


9684 








-0.91 ± 0.09 


0.95 


9845 








-1.12 ±0.06 


0.95 


10050 


-1.80 ± 0.18 


0.98 


-4.03 ± 0.29 0.95 






10247 








-1.27 ±0.07 


0.95 


10306 








-1.27 ±0.09 


0.94 


10323 








-0.98 ± 0.05 


0.95 


10488 


-1.60 ±0.16 


0.94 


-3.64 ±0.19 0.94 






10570 








-0.83 ± 0.07 


0.95 


MEAN 


-1.90 




-4.32 


-1.01 




o~- 


0.35 




0.86 


0.36 





part. Both the peak and total energy distribution functions follow 
single power laws with indices -1.63 + 0.15 and -1.45 + 0.13, 
respectively. Figure 5 depicts the distribution functions of dura- 
tion (Fig. 5a), peak energy (Fig. 5b) and total energy (Fig. 5c) for 
AR9415. Here the duration distribution follows a power law with 
an exponential rollover. The power law index is - 1 .42+0. 1 8. The 
peak and total energy distribution functions follow again single 
power laws with indices -1.84 + 0.18 and -1.50 + 0.13, respec- 
tively. 

From the above, flare-duration distributions appear to be 
best-fitted either by power laws or by power laws with expo- 
nential rollovers. By comparing the chi-square values and re- 
spective probabilities for single power laws, double power laws 
and power laws with exponential rollover respectively, it is con- 
cluded that 3 of our ARs follow double power laws (AR9635, 
AA10050, AA10488), while the rest are best fitted by power 
laws with exponential rollovers. These indices are summarized 
in Table 1, along with the respective probabilities. The values 
shown in Table 1 refer to the fitting achieved against the en- 
tire distribution function in all cases (all bins included). Single 
power laws fail to describe the model duration distribution func- 
tions in all cases, whereas exponential functions only fit the tail 
of the generated model curves. In cases where a double power 
law is the best fit (AR9635, Atf 10050, AA10488), the mean in- 
dex for the flat power law is - 1 .90, whereas the mean index for 
the steep power law is -4.32. When power laws with exponential 
rollover are best fitting (remaining ARs), then the mean value for 
the power law index is -1.01. Standard deviations (cr 2 ) to these 
mean values are given in the last row of this table. 

Although single power laws are not the optimum functions 
to fit the modeled flare duration, they are undoubtedly the best- 
fitting theoretical functions for the peak energy and the total flare 
energy. As shown in Table 2 for the peak flare energy and in 
Table 3 for the total flare energy, the average value for E pea k is 
-1.80, whereas the average index value for E, ota i is -1.57. The 
standard deviation (cr 2 ) of these mean values is given in the last 
row of these tables. 

Figure 6 illustrates the magnetic energy released E re [ for a 
specific period of 10000 timesteps after SOC has been reached 
for A/?10570. As the added driver increments SB assume small 
and random values, the waiting time from one flaring event to 
another varies. Figures 7 and 8 show a 3d representation of the 
emerging magnetic discontinuities during a large and a smaller 



Table 2. Power-law indices and respective chi-square probabili- 
ties derived by fitting a single power law to the peak flare energy 
for the 1 1 ARs comprising our sample. 



FLARE PEAK ENERGY (model) 
Single Power Law fit 



AR 


PL Index 


Probability 


9415 


-1.84 + 0.18 


0.95 


9635 


-2.62 ±0.17 


0.97 


9661 


-1.42 ±0.15 


0.98 


9684 


-1.70±0.17 


0.97 


9845 


-1.85 ±0.12 


0.95 


10050 


-1.63 ±0.15 


0.95 


10247 


-2.15 ±0.12 


0.98 


10306 


-1.61 ±0.16 


0.97 


10323 


-1.72 ±0.17 


0.97 


10488 


-1.59 ±0.14 


0.95 


10570 


-1.63 ±0.15 


0.98 


MEAN 


-1.80 




T 

cr 


0.33 





avalanche, respectively, simulated for AT? 10247 after SOC has 
been reached. In the former case, the avalanche during its early 
stages generates 140 discontinuities (Fig. 7a), evolves further 
(Fig. 7b) with 281 discontinuities, peaks (Fig. 7c) with 425 
discontinuities, and decays (Fig. 7d, 7e, 7f) with 184, 51, and 
18 discontinuities, respectively. The total event duration is 341 
steps. The total duration of the smaller event (Fig. 8) is 90 steps. 
The event during its early stages generates 6 discontinuities (Fig. 
8a), peaks (Fig. 8b) with 15 discontinuities, and decays (Fig. 8c, 
8d) with 10 and 6 discontinuities, respectively. 

Finally, it is interesting to investigate whether our 
model consistently reproduces the distribution func- 
tions of the flaring events actually observed in the 
ARs in our sample. For our comparison we used 
the solar X-ray flare catalog from the GOES satellite 
(http : //www . ngdc . noaa . gov/stp/SOLAR/f tpsolar flares . html 
item 3). The flaring events recorded in this database lie in the 
class range B - X and are summarized in Table 4 for each 
AR. However, to construct the distribution functions of flare 
parameters we need sufficient statistics, reflected in large flare 
numbers. A single AR, regardless of flare productivity, is 
unlikely to provide these numbers. For this reason and for the 
sake of comparison we have merged all observed flares in all 
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Fig. 4. Distribution functions for the event duration (Fig. 4a), the peak energy E pea i (Fig. 4b) and the total energy E tota i (Fig. 4c) for 
ARIQQ5Q. The energies in b) and c) are calculated in physical units (ergs). 



studied ARs into a single flare sequence with a total of 154 
events (sum of all flares in Table 4). Table 5 shows the statistical 
results of our analysis for flare durations. As flare duration we 
define the observed onset-to-end elapsed time. The best fitting 
function is not easily discernible in this case, as all candidate 
functions (single power law, double power law and power law 
with exponential rollover) fit the observational data fairly well. 



Figure 9 depicts the fit between the observed flare durations 
against a double power law (Fig. 9a) and a power law with 
exponential rollover (Fig. 9b). In the former case the calculated 
index is -1.67 + 0.09 for the smooth part and -3.37 + 0.25 for 
the steep part, whereas in the latter case the power law index 
yields -1.28 ± 0.11. It is apparent that the dynamical range of 
the power-law in Fig. 9a is very limited, but it is shown here 
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Fig. 5. Same as Fig. 4 for AR9415. 



for comparison purposes. In order to achieve this comparison, 
we merge the model results of the separate runs per AR into 
one common database. Figure 10 depicts the fit between the 
merged model flare durations for all ARs in our sample against a 
double power law (figure 10a) and a power law with exponential 
rollover (Fig. 10b). In the former case the calculated index 
is -1.78 ± 0.27 for the smooth part and -3.91 ± 0.42 for the 
steep part, whereas in the latter case the power law index 
yields -1.13 ± 0.12. By comparing the results depicted in Fig. 



9 (observational data) and Fig. 10 (merged model data), we 
conclude that the power law indices for the observational data 
are close to our model's values for the double power law and 
power law with exponential rollover fits. This is not the case for 
the single power law fitting, which yields an index of - 1 .70 for 
the GOES data. The best agreement between the observed and 
the simulated flares is, therefore, achieved when the attempted 
fit is not a single power law, but either a double power law or 
a power law with an exponential rollover. Table 6 is similar 
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Fig. 6. Time series of the total magnetic energy released E re i (in 1 .033xl0 25 ergs) from the simulated A/?10570 after the SOC state 
has been reached. The time series shown consists of 10000 timesteps for better detail. 



Table 3. Power-law indices and respective chi-square probabili- 
ties derived by fitting a single power law to the total flare energy 
for the 1 1 ARs comprising our sample. 

FLARE TOTAL ENERGY (model) 



Single Power Law fit 



AR 


PL Index 


Probability 


9415 


-1.50 ± 0.13 


0.95 


9635 


-2.22 ±0.19 


0.98 


9661 


-1.27 ±0.05 


0.99 


9684 


-1.43 ±0.07 


0.99 


9845 


-1.69 ±0.17 


0.95 


10050 


-1.45 ±0.13 


0.95 


10247 


-1.89 ±0.17 


0.98 


10306 


-1.23 ±0.08 


0.99 


10323 


-1.45 ±0.16 


0.98 


10488 


-1.54 ±0.13 


0.95 


10570 


-1.45 ±0.08 


0.99 


MEAN 


-1.56 




cr 1 


0.28 





Table 4. GOES X-ray data for the number and class of observed 
flares in the ARs used in our simulations. 



AR B-Class C-Class M-Class X-Class Total 



9415 


03 


16 


06 


05 


30 


9635 


00 


02 


00 


00 


02 


9661 


00 


16 


01 


02 


19 


9684 


00 


08 


01 


01 


10 


9845 


01 


04 


00 


00 


05 


10050 


00 


16 


00 


00 


16 


10247 


00 


01 


00 


00 


01 


10306 


06 


02 


00 


00 


08 


10323 


00 


05 


00 


00 


05 


10488 


00 


17 


07 


02 


26 


10570 


17 


14 


01 


00 


32 



to Table 5, summarizing the indices resulting from the merged 
model data. 

Although our findings show good alignment both with pre- 
vious models and observations, it is crucial to crosscheck the 
physical soundness of our algorithm. As mentioned in Sect. 3.4, 
loading rule (10) does not by itself guarantee that the mag- 
netic field remains divergence-free during the entire simulation. 



WNDB is therefore determined in order to monitor the mag- 
netic field divergence throughout the loading and redistribution 
process. Figure 1 1 presents the evolution of WNDB during the 
3 x 10 5 timesteps of our simulation for AR1024-7. In the begin- 
ning WNDB is close to zero, as our initial condition is the ex- 
trapolated NLFF (and therefore approximately divergence-free) 
magnetic field. As time elapses, V • B starts deviating from zero, 
but WNDB remains under 0.20 during the entire simulation. This 
holds for all ARs in our sample. Therefore, our model retains 
the magnetic field approximately divergence-free throughout the 
simulation. 

Furthermore, it is worth investigating whether the use of 
alternative threshold definitions incurs any qualitative changes 
in the presented results. As an example, we apply the sec- 
ond threshold definition of Sect. 3.2 to ARIQ2A-1 . In this case 
G cr = (?ai W (l - s) =s 30G. Figure 12 shows that even with 
this threshold definition, G av increases gradually until timestep 
1.3xl0 5 , after which the SOC state is reached. G av stabilizes 
around approximately 29.95, which is lower than the critical 
threshold G cr = 30. The statistical properties of the generated 
distribution functions remain unchanged. This is also valid when 
switching from a 32 X 32 X 32 grid towards larger volumes (, e.g. 
a 64 x 64 x 64 grid). 

5. Discussion and conclusions 

This study simulates the flaring activity of 1 1 solar ARs in terms 
of a refined CA model. The modules comprising this integrated 
flare model are summarized below: 

1. We extrapolate the magnetic field from the photospheric 
boundary of 11 IVM magnetograms resampled on a 32 X 
32 grid, through a nonlinear force-free optimization algo- 
rithm with preprocessing at the photospheric level (module 
"EXTRA", refer to Sect. 3.1 for details). 

2. We identify the unstable locations which will dissipate mag- 
netic energy in our grid when the approximated magnetic 
field Laplacian G av (r) at site r exceeds a specific threshold 
G cr (module "DISCO", refer to Sect. 3.2 for details). 

3. In case magnetic discontinuities are identified (either directly 
after the initialization or after each loading), the magnetic en- 
ergy is redistributed such that the instabilities are completely 
relaxed (module "RELAX" refer to Sect. 3.3 for details). 
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Fig. 7. 3d representation of the emerging magnetic discontinuities during an avalanche in A/?10247. The total duration of this event 
is 341 steps. During the early stages, the avalanche generates numerous discontinuities (140 in a), evolves with 281 discontinuities 
(b), peaks with 425 discontinuities (c), and decays with 184 discontinuities (d), 51 discontinuities (e), and 18 discontinuities (f). 



4. Further loading within our system is allowed, when a pre- 
viously triggered avalanche has completely decayed. In this 
case, LOAD adds a random magnetic field increment 6B(r) 
at a random site r within our grid according to the rules de- 
scribed in Sect. 3.4. 



The algorithm is allowed to run for 3xl0 5 timesteps, which 
is sufficient for all simulated ARs to both reach the SOC state 
and provide sufficient event statistics after the SOC state has 
been reached. The enhancements of our flare simulation model 
in comparison with previous SOC models of solar flares are fol- 
lowing: 



Dimitropoulou et al.: Simulating flares in active regions driven by observed magnetograms 



11 



a) 



C) 





b) 



d) 





Fig. 8. 3d representation of the emerging magnetic discontinuities during an avalanche in A/?10247. The total duration of this event 
is 90 steps. During the early stages, the event generates a small number or discontinuities (6 in a), peaks with 15 discontinuities (b), 
and decays with 10 discontinuities (c) and 6 discontinuities (d). 

Table 5. Power law indices and respective chi-square probabilities derived by fitting several functions to the flare durations derived 
from the merged GOES observational data for the 1 1 ARs comprising our sample. 



FLARE DURATION (data) 



Single Power Law 



Double Power Law 



Power Law w Exponential Rollover 



Flat Power Law 



Steep Power Law 



PL Index Probability PL Index Probability PL Index Probability PL Index 



Probability 



-1.70 ±0.12 



0.98 



-1.67 ±0.09 



0.98 



-3.37 ±0.25 



0.95 



-1.28 ±0.11 



0.96 



Table 6. Power law indices and respective chi-square probabilities derived by fitting several functions to the flare durations derived 
from the merged model data for the 1 1 ARs comprising our sample. 



FLARE DURATION (merged model data) 



Single Power Law 



Double Power Law 



Power Law w Exponential Rollover 



Flat Power Law 



Steep Power Law 



PL Index Probability PL Index Probability PL Index Probability PL Index 



Probability 



-2.79 ± 0.22 



0.97 



-1.78 ±0.27 



0.98 



-3.91 ±0.42 



0.94 



-1.13 ± 0.12 



0.96 



The initial boundary conditions are not arbitrary, but stem 
from real solar magnetograms. An NLFF field extrapolation 
is used to reconstruct the initial magnetic configuration gen- 
erated from the observed 1 1 ARs, retaining to the best pos- 
sible extent physical requirements such as the minimization 
of the Lorentz force and the magnetic field divergence. 



- Given that the simulation commences from observed mag- 
netograms, it is now possible for our CA model to remove 
the restriction of arbitrary energy units (see e.g. the remarks 
within Georgo uiis et al. (2001)| l. This gives us the opportu- 
nity to directly compare the model with the observed energy 
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Fig. 9. Observed distributions of GOES flare durations for all ARs in our sample. Fit is attempted using a double power law (a) and 
a power law with an exponential rollover (b). The double power law fit yields an index equal to -1.67 ± 0.09 for the flat part and 
-3.37 ± 0.25 for the steep part, whereas the fit with the power law and the exponential rollover yields a scaling index - 1 .28 ± 0. 1 1 . 
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Fig. 10. Distributions of simulated flare durations for all ARs in our sample. Fit is attempted using a double power law (a) and a 
power law with an exponential rollover (b). The double power law fit yields an index equal to —1.78 ± 0.27 for the flat part and 
-3.91 + 0.42 for the steep part, whereas the fit with the power law and the exponential rollover yields a scaling index -1.13 ± 0.12. 




Fig. 11. Evolution of WNDB during the 300000 timesteps of our simulation for AA10247. 
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Fig. 12. Diagram of the average G av value over the grid for 3xl0 5 timesteps for AR102A1 when the threshold definition is G cr = 
G aVm Jl - s) - 30G. 



content per flare, thus leaving time as the only arbitrary quan- 
tity in our simulation. 

Our model follows to a significant degree the principles of 
|Lu & Hamilton (1991)| The rules obeyed during both the 
magnetic energy redistribution and the further driving of the 
system are designed in such a way that the magnetic field di- 
vergence is within tolerated limits. This has not been the case 
in the early CA models (Vlahos 1995, Georgoulis & Vlahos 
1996, Georgoulis & Vlahos 1998) and has only been touched 
in advanced C A approaches through the use of the vector po- 
tential A instead of the magnetic field B in combination with 
an improved way of calculating the derivatives (Isliker et al. 
2000, Islikeretal. 2001). 

The driving mechanism attempts to mimic not only pho- 
tospheric convection as proposed by Parker (1988, 1989, 
1993), but also coronal evolution, such as turbulence and 
current sheet interaction. In this sense, locations through- 
out the simulation box are randomly chosen to be perturbed. 
Turbulence, via either localized Alfven waves or larger-scale 
turbulent flows (Einaud fet al. 19961 Rappazzo et al. 2008 ) 
leads to current-sheet interaction that, depending on the 
local magnetic conditions, may trigger an avalanche ob- 
served as a flare. Naturally, though, due to the larger ac- 
cumulation of magnetic free energy close to the photo- 
sphere ( Regnie r"& Priest 2001) photospheric convection and 
systematic photospheric motions (e.g. shear) should be the 
drivers of most coronal instabilities. Our driving mecha- 
nism should be further revised to account for systematic 
photospheric flows. This future step is important because 
(1) it has been argued that the distribution and energy con- 
tent of magnetic discontinuities in a given photospheric 
boundary can explain the statistical properties of flares 
( Vlahos & Georgoulis 2004) and (2) investigating possible 
correlations between the photospheric driver and the cor- 
responding coronal active region reveals the strong nonlin- 
earity of active-region magnetic configurations that hinders 
correlations between photospheric and coronal structures 
(Dimitropoulou et al. 2009). The latter patterns, however, 
have a crucial impact on the expected dynamical activity of 
the system, namely, the magnetic energy release and the sub- 
sequent particle acceleration processes (Vla hos et al. 2004b . 
The derived results can be directly compared with flare ob- 
servations, due to the fact that the simulation uses extrap- 
olated fields from observed vector magnetograms as initial 



conditions. At this point, we once again stress that the X-ray 
flares recorded by GOES for each AR do not comprise a sta- 
tistically reliable sample. Therefore, in order to make such a 
comparison possible, we merged the GOES flare data of all 
1 1 ARs in our sample into one database, comprising of 154 
flares. 

Our results show that under the imposed driving and re- 
distribution rules, all examined ARs reach the SOC state. The 
retrieved distribution functions for event duration are best de- 
scribed by either double power laws or power laws with expo- 
nential rollover, although single power laws are also applicable 
for the merged data. The peak energy and total energy follow 
clearly single power laws. The power law indices for durations 
and energies as presented in Tables 1,2 and 3 lie in the well- 
known ranges documented consistently in numerous past stud- 
ies, including[Georgouli s et al. (2001)| In this study, Georgoulis 
et al. compare their SOC model with data from the Danish 
Wide Angle Telescope for Cosmic Hard X-rays (WATCH) col- 
lected during maximum of the solar cycle 21. Figure 1 in the 
cited work shows that the peak and total energy of the observed 
flares follow single power law distribution functions with in- 
dices -1.59 and -1.39 correspondingly, whereas the flare du- 
ration distribution function is considered to either follow dou- 
ble power law (with index -1.15 for the flat and -2.25 for the 
steep part) or power law with exponential rollover (with power 
law index -1.09). These results are in agreement with our find- 
ings. Although our model generates flare duration distribution 
functions with indices in alignment with the ones presented in 
|Georgoulis & Vlahos (1998)| we did not attempt here t o repro- 
duce two key findings of Georgoulis & Vlahos (1998) namely 
the variability of the scaling indices as a result of the driver's 
variability as well as the two distinct event populations. In the 
cited study the peak and total energy distribution functions fol- 
low double power laws. The steeper part of them corresponds to 
the signature of a "soft" flare population (nanoflares), whereas 
the flatter part is attributed to microflares and flares. 

Although this work overcomes major drawbacks of many 
previous CA models, such as retaining the value of the mag- 
netic field divergence close to zero throughout the simulation, 
there are still some points that can lead to discrepancies. First 
and foremost, the determination of the threshold value G cr can 
slightly influence the exponents of the retrieved power laws, al- 
though it cannot cause any qualitative change to their appear- 
ance, namely the known flare statistical properties will always 
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follow power law distributions, independent of the threshold 
value imposed. The histogram method presented in Sect. 3.2 
eliminates to an extent the arbitrary selection of G cr . We have 
also investigated whether the rebinning of our grid to the size of 
32x32x32 influences our results in comparison with larger grids 
(e.g. 64 x 64 x 64) and we found that the differences in the power 
law indices lie within the inferred uncertainties. Finally, as far 
as the comparison with the observational data is concerned, we 
have already stressed that this is a preliminary attempt given that 
the number of GOES X-ray flares across all investigated ARs 
does not produce sufficient statistics. 

The discussion regarding the validity of the CA models when 
it comes to the simulation of physical processes in complex sys- 
tems is a long-running one. As discussed by Islik er et al. (1 998), 
the essence of CA modeling is to describe complex systems, 
which comprise a large number of interacting subsystems, as- 
suming that the global dynamics described statistically are not 
sensitive to the fine structure of the elementary processes. More 
strict approaches such as MHD, on the other hand, are based on 
the precise description of the elementary processes through de- 
tailed differential equations. Both approaches have been shown 
to exhibit drawbacks and advantages. The CA approach does not 
provide any insight into the local processes or over short time in- 
tervals, but it reproduces the global statistics. MHD reveals de- 
tails about the local processes, but coupling them to a global de- 
scription is a formidable task. In this sense, the two approaches 
are complementary and there have been indeed various attempts 
to either combine them (e.g. Longope & Noonan 2000| l, or in- 
terpret CA models as discretized MHD equations (Isliker et al. 
1998,Vassiliadis et al. 1998). Even more extended CA models, 
like the X-CA model described by Isliker et al. (2001), have 
achieved consistency with MHD to a greater extent. Our CA 
model will opt to incorporate and utilize meaningful modeling 
developments into a more concrete, "integrated" flare model. 
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